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Abstract: 

Previous studies have drawn attention to substantial hydrological changes taking place in mountainous watersheds where hydrology is 
dominated by eryospheric processes. Modelling is an important tool for understanding these changes but is particularly challenging in 
mountainous terrain owing to scarcity of ground observations and uncerliiinty of model parameters across space and lime. This study 
utilizes a Markov Chain Monle Carlo data assimilation approach to examine and evaluate the performance of a conceptual, degree-day 
snowmelt runoff model applied in the Tamor River basin in the eastern Nepalese Himalaya. The snowmelt runoff model is calibrated 
using daily streamflow from 2(X)2 to 2006 with fairly high accuracy (average Nash Sutcliffe metric -0.84, annual volume bias < 3%). 
The Markov Chain Monte Carlo approach constrains the parameters to which the model is most sensitive (e.g, lapse rale and recession 
coefficient) and miiximizes model fit and performance. Model simulated stieamllow using an interpolated precipitation data set 
decreases the fractional contribution from rain tail compared with simulations using observed station precipitation. The average 
snowmelt contribution to total runoff in the Tamor River basin for the 2(X)2 21X16 period is estimated to be 29.7 ±2.9% (which includes 
4.2 ±0.9% from snowfall that promptly metis), whereas 70.3 ±2.6% is attributed to eonlribu lions from rainfall. On average, the 
elevation zone in the 4000 5500m range contributes the most to basin runoff , averaging 56.9 ±3.6% of all snowmelt input and 
28.9 ± 1.1% of all rainfall input lo runoff. Model simulated stieamllow using an interpolated precipitation data set decreases the 
fractional contribution from rainfall versus snowmelt compared with simulations using observed station precipitation. Model 
experiments indicate that the hydrograph itself does not constrain estimates of snowmelt versus rainfall contributions to total outflow 
bul that this derives from the degree-day melting model. Lastly, we demonstrate that the data assimilation approach is useful for 
quantifying and reducing uncertainly related to model parameters and thus provides uncertainty bounds on snowmelt and rainfall 
contributions in such mountainous watersheds. Copyright © 2013 John Wiley & Sons, Ltd. 

Supporting information may be found in the online version of this article. 
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INTRODUCTION 

Growing evidence indicates that the low-latitude mountain- 
ous cryosphere is undergoing change owing to recent 
warming (IPCC, 2007; Yao et al . , 2012). Geographic areas 
where Lhe water cycle is dominated by snowmell hydrology 
are expected lo be more susceptible lo climate change, as il 
affects the seasonality of runoff (Barnett et al. , 2005; Adam 
et al . , 2009). The Hindu Kush-Himalaya (HKH), regarded 
as lhe H water towers 1 of Asia, is one such critical region as it 
nourishes several of the major Asian watersheds through 
perennial rivers such as the Ganges, Indus and Brahmaputra 
(Messerli et al., 2004). Most recent research analyzing 
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tropospheric temperatures from the longest available record 
of microwave satellite measurements reveals widespread 
pre-monsoonal and mean annual warming over the HKH 
region at a rate of 0.21 ±0.08 °C/decade (Gautam et al., 
2010). Regional climate projections by lhe Intergovern- 
mental Panel on Climate Change (IPCC, 2007) indicate 
Central Asia may warm by a median temperature of 3.7 °C 
by the end of the 21st eenluiy, with warming over high 
altitudes across the HKH (Panday et al ., in review). There is 
an increasingly large body of evidence of a glacio- 
hydrological response along the east-west transect of the 
HKH corresponding to the these climatic changes, with 
glaciers in the Eastern Himalaya exhibiting retreat and 
negative mass balance and glaciers in the Karakoram and 
northwestern Himalaya exhibiting a positive mass balance 
over the last few decades (Bolch et aL , 2012; Yao et al., 
2012). Such climate-driven responses of mountainous river 
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hydrology may pose significant challenges for this region. 
Despite its regional importance and vulnerability to climate 
change, there is uncertainty associated with the physical 
processes that control runoff generation at high altitudes 
(Boleh et al. , 2012; PelliceioLti et al. , 2012) and subse- 
quently with the rates and magnitude of climate change 
impacts on snow cover and snowmelt hydrology. 

Water resource management and the evaluation of 
impacts of climatic change require quantification of 
stream flow variability. Hydrological models provide a 
framework to investigate these relationships and subse- 
quently use a scenario approach to evaluate projections of 
hydrological variables (Maurer et al ., 2009). The assessment 
of hydrological impacts of climate change is particularly 
challenging in mountainous environments as Lhey constitute 
complex hydrological systems owing to extreme heteroge- 
neity in vegetation, soils, topography and spatially and 
temporally varying snow cover and snowmelt patterns (Gurtz 
et a/., 2005; Panday et al ., 2011). Continued efforts in 
hydrological modelling of snow and glacier melt in the HKH 
region mostly through conceptually lumped, semi -distributed 
(rather than physically based models) have provided baseline 
patterns of spaliotemporal patterns of precipitation and runoff 
at very large scales (Bookhagen and Burbank, 2010; 
Immerzecl et al ., 2010). However, poor accessibility and 
scarcity of adequate network of hydrometeorological stations 
at high altitude regions still remain a major impediment to 
runoff modelling. As a result, only a few studies have 
explored snowmelt runoff models al the catchment or river- 
basin scale in this region (Braun et al. y 1993 ; Shreslha et al., 
2004; Akhlar et al. , 2008; Tahir et al., 2011). 

Parameterization is usually challenging owing to scarcity 
of data in mountainous regions such as the HKH, and a 
sound understanding of the snowmelt hydrology needs to be 
emphasized prior to assessment of water resources in a 
changing climate. Several catchment-scale and basin-scale 
modelling studies have addressed the glacio-hydrological 
processes in the Himalayan region (Konz et al. , 2007; 
Akhlar et al ., 2008; Butt and Bilal, 201 1; Immerzeel et al., 
2011; Tahir et al., 2011). A majority of these studies have 
adopted parameter values from the literature without 
rigorous calibration and without quantifying how parameter 
uncertainty influences modelled results. Transfer of param- 
eters across space and time in hydrological modelling can be 
problematic, particularly in conceptual models of snowmelt 
when parameters either do not have a physical meaning or 
Lire not easily measurable (Sun et al., 2012). In order to 
address the issues of parameterization and model uncertain- 
ly, data assimilation methods are typically applied in 
hydrology to integrate different datasets, combine model 
outputs with observations, update hydrological models, 
quantify model uncertainties and sensitivity and/or imp rove 
model accuracy (Beven and Freer, 2001; Clark etal. , 2006). 
Data assimilation can be viewed as a method to use 


observations to constrain models with the goal of providing 
better understanding of the underlying processes 
(McLaughlin et al ., 2005). The Bayesian approach in 
hydrological applications allows to systematically incorpo- 
rate prior information on model parameters in the data 
assimilation framework (Zobilz etai, 2011) and to compute 
probability distributions of modelled outputs and character- 
ize uncertainty in the model, usually through Monte Carlo 
simulation (Mishra, 2009). This provides a means of 
objectively and quantitatively diagnosing the parameters 
and processes to which a model is most sensitive and from 
which one can learn. Few studies have emphasized the need 
for parameterization in the Himalayan region and have 
performed extensive calibration (Konz et al ., 2007; 
Immerzeel et al ., 2011; Pellicciotti et al ., 2012) and 
validation with additional data sets such as glacier mass 
balances and remotely sensed snow cover. Pellicciotti et al. 
(2012) also emphasized that an assessment of the internal 
consistency of the model is important owing to the problem 
of multiplicity of parameter sets or equifinalily. 

In addition to parameter uncertainty, sensitivity of 
hydrological models to input data such as precipitation 
also needs to be evaluated particularly in complex 
topography and high elevations. Hydro logi cal models in 
mountainous regions with scarce observations arc usually 
forced with precipitation data located at lower altitudes, 
which can introduce added uncertainly. Satellite-derived 
precipitation estimates and interpolated rain-gauge obser- 
vations have been used in other studies as inputs to a 
snowmelt runoff model (SRM) instead of gauged precip- 
itation. Bookhagen and Burbank (2010) used improved 
Tropical Rainfall Measuring Mission rainfall estimates to 
characterize the spatiotemporal distribution of rainfall, 
snowfall and evapotranspi ration across the Himalayan 
region, whereas Tahir et al. (201 1) successfully used 
interpolated daily precipitation data to simulate snowmelt 
hydrology in the Hunza River basin. 

A thorough assessment of hydrologic processes, 
patterns and trends requires minimization of uncertainties 
in modelling and available input data. The objective of 
this study is to examine spatial and temporal variations in 
snowmelt contributions to runoff in the Hindu Kush- 
Himalayan region and Tamor River basin by using a data 
assimilation approach that both combines available input 
data with a degree-day based SRM and estimates model 
parameters by using a Markov Chain Monte Carlo 
(MCMC) technique. Additionally, we investigate and 
evaluate the performances of the SRM when driven by 
observed precipitation as well as interpolated gridded 
precipitation data. We demonstrate how a data assimila- 
tion approach can be coupled with a SRM to investigate 
the issues of model parameter sensitivity and uncertainly 
and improve performance in a data-scarce basin in the 
Eastern Himalaya. 
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METHODS 

Study area 

This geographic domain of Lhis study is the upstream 
portion of the Tamor River basin in the Eastern Himalaya 
of Nepal (Figure 1). The Tamor basin covers an area of 
4274 km 2 , with an elevation range of 422 to 8505 m 
(Figure 2). The Tamor River, which is one of the major 
tributaries of the Sapta Koshi River, flows between high 
peaks such as Kanchenjunga (8586 m) and Everest 
(8848 m) on gorges reaching depths of 2000 m. Precip- 
itation for the region is highly seasonal and concentrated 
during the summer monsoon months (June to September) 
(Figure 3). The upper reaches of the basin are snow 
covered and glaciated, contributing melt water to 
stream flow. 

Hydrometeorological data 

Hydrometeorological measurements for the Tamor 
RiveF basin in the Eastern Himalaya were acquired from 
the Department of Hydrology and Meteorology, Nepal. 
Daily temperature data at Taplejung station (27,35 °N, 
87.67°E at 1732m elevation) and daily precipitation data 
measured at Lungthung station (27.55°N f 87.78°E at 
1 780 m elevation) in the Taplejung district were used for 
the hydrological modelling (Figure 1). Discharge data 
measured daily from 2002-2006 for the Tamor River at 
Majhilar station (27. 1 5°N, 87.7 1 °E at 533 m elevation) were 
used for calibration. The Tamor basin was delineated at the 
Majhilar station using the 90m digital elevation data 



acquired from the Consultative Group on International 
Agricultural Research - Consortium for Spatial Information. 
The basin was divided into four zones, and the zonal mean 
hypsometric elevation (Figure 2) was used as the elevation 
to which base station temperatures were extrapolated for the 
calculation of zonal degree days (Table l). 

Moderate Resolution Imaging Spectroradiometer (MOD IS) 
snow cover 

Satellite-derived snow cover is the most efficient and 
readily available input for snowmelt runoff modelling. 
Snow cover was derived using the MODIS/Terra Snow 



Rivers 

Meteorological stations 
® Majhitar station 


Elevation 


1H 01500 

HI 1500-2100 
I 1 2100-2800 
| 2800-3400 
k— I 3400-3900 
EMI 3900-4500 
^B 4500-5100 

^B 5100-5600 

HI 5600-6400 
! 6400-8X00 


Figure L Location and topography of the Tamor River basin in the eastern Nepalese Himalaya 
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Figure 3. Average monthly slreamflow (m 3 /s) and total monthly 
precipitation (mm) from 1996 to 2006 at Majhitar station 


Asian Precipitation □ Highly Resolved Observational 
Data Integration towards Evaluation (APHRODITE) 
precipitation data 



used to drive a snowmelt runoff model, as valley stations 
may not be representative for basin precipitation. The 
APHRODITE of the water resources project has compiled a 
daily gridded precipitation dataset covering 1951-2007 
period by analyzing rain-gauge observation data across Asia 
(Yatagai et al., 2009). The APHRODITE data used in this 
research (APHRO_MA_Vt003Rt) have a 0.25x0.25° 
gridded spatial resol ution, which was interpolated from 
meteorological stations throughout this region. Whereas 
satellite data consistently underestimate or overestimate the 
amount of precipitation, particularly in the monsoon season 
(Duncan and Biggs, 2012), APHRODITE data have been 
shown to provide the best precipitation estimates in most 


Table I. Elevation zone, elevation range, mean elevaLion and 
zonal area for the four zones of Lhe Tamor River basin 


ElevaLion zone 

Elevation 
range (m) 

Mean elevaLion 
of zone (m) 

Zonal 
area (km 2 ) 

1 

533-2500 m 

1685.9 

1355.5 

2 

2500 40(X)m 

3208.1 

1111.3 

3 

4000 5500 m 

4764.3 

1297.9 

4 

>5500 

5988.8 

509.2 

Total basin 

— 


4273.8 


Cover 8-Day L3 Global 500 m Grid (MOD10A2; Hall etal . , 
2006) and updated every eight days using a snow depletion 
curve. The MODIS snow cover is based on a snow mapping 
algorithm that uses at-s ate Hite reflectance in MODIS bands 
4 (0.545-0.565 in) and 6 (1 .628-1 .652 m) to compute the 
Normalized Snow Difference Index (Hail et al 2002). A 
pixel will be mapped as snow if the Normalized Snow 
Difference Index U 0.4 and MODIS band 2 (0.841- 
0.876 m) >11%; however, if the MODIS band 4 
reflectance is < 10%, the pixel will not be mapped as snow 
even if the other criteria Lire met. The MODIS product 
distributed through Lhe National Snow and Ice Data CenLre 
was acquired for the scene h25v06, covering the study 
region from 2002 to 2006. The MODIS images were 
reprojected to Universal Transverse Mercator zone 45 
projection system, and the Tamor River basin area was then 
extracted to assess the fraction of snow cover for different 
altitudinal zones. Average snow cover was estimated by 
interpolating linearly between the previous and next 
available cloud free images for days with cloud cover 
following Tahir et al. (201 1). The conventional depletion 
curves were derived for each elevation zone using the 8 -day 
snow cover data, which provide the input to the SRM as the 
decimal portion of snow-covered area for each day. 


set is limited (Andermann etal., 2011). We used the basin- 
averaged APHRODITE precipitation as input to the SRM 
assuming it to be constant across elevation zones. The model 
was run using APHRODITE preeipitadon data to determine 
the reliability of the dataset for snowmelt mnoff modelling 
in the Tamor River basin. 

Snowmelt runoff modelling 

Modelling of stream discharge was based on the SRM 
developed by Martinec et al. (2008). SRM is a conceptual, 
detenmnistic, degree-day hydrological model that simulates 
and forecasts daily runoff resulting from snowmelt and 
precipitation from hydrometeorological input data. The model 
was designed for streamfiow modelling in basins where 
snowmelt is a major contributor, and more recently, it has 
been applied for evaluation of impacts of climate change on 
seasonal snow cover and runoff (Immerzeel et al ., 2010). 
SRM has been applied to basins ranging from 0.76 to 
917 444 km 2 and across elevation ranges from 0 to 8840 m in 
over 100 basins in 29 different countries (Martinec et al ., 
2008). This model was chosen in this study to assess 
snowmelt hydrology owing to its minimal data requirements 
and transferability to a variety of geographic regions. On the 
basis of the input of daily temperature, precipitation and snow- 
covered area values, SRM computes the daily streamilow as 

Gniii □ □ E LlT m uS m (1) 

UCr m P m 1A Z 1 (XXX) / 86400 “□ □ A„ ul □ 

— Qn L 1 1 

*„oi □ (2) 

where Q (m 3 /s) is discharge al day n + 1 , C is lhe runoff 
coefficient lo snowmelt (Cv n ) and to rain (Cr„) for each zone 
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(subscript z), a (cm/°C/day) is the degree-day factor, 
T + LfT (°C) arc the degree days, P (cm) is the precipitation 
per day, A (km 2 ) is the area of the basin, S is the fractional 
snow cover and k is the discharge recession coefficient. The 
daily average stream flow on day n + 1 is computed by 
addition of snowmelt and precipitation that contributes to 
runoff with the discharge on the preceding day. The 
snowmelt of the preceding day is a product of the degree- 
day factor (a), zonal degree days (T + _T) and the snow- 
covered area percentage (S). The snowmelt runoff coefficient 
( Cs ) and the zonal area (A) arc further multiplied with the 
previous product to compute the percentage that contributes 
to runoff. Similarly, the product of the rainfall runoff 
coefficient (Cr) and the zonal area (A) determines the 
precipitation contribution to runoff. The recession coefficient 
(£) is also an important parameter as it describes the structure 
of the falling limb of the hydrograph and correspondingly 
detemiines the proportion of daily melt water production that 
appears immediately in the runoff. This is obtained with 
analysis of historical discharge data based on Equation 2, 
where parameters x and y are derived by analyzing discharge 
on a given day as a function of the value on the previous day. 

A critical temperature is used to determine whether a 
precipitation event is rain or snow over different parts of the 
basin. When precipitation is determined to be snow, its 
contribution to runoff depends cm whether it falls on a snow- 
covered or previously snow- free area. Snowfall on snow- 
covered areas is assumed to become part of the seasonal 
snowpack, its depth is not tracked further, and its 
contribution to runoff is determined by the snow depletion 
curve. Daily snowmelt depth from snow-covered areas is 
determined on the basis of degree days and a factor that 
converts this to mass release, qualitatively representing the 
maximum potential melt given the energy input. In contrast, 
snow failing as precipitation over snow-free areas is 
accumulated until the next day warm enough to produce 
melting at which time it is depleted on the basis of the 
degree -days approach. During the winter half year 
(October— March), the model cancels any existing storage 
of preceding temporal snowfalls when snow cover is 
maximum such that snow on previously snow- free areas 
becomes part of the seasonal cover. When precipitation is 
considered to be rain during the cold season, rainfall runoff 
contributes to total runoff only from the snow- free area, and 
the rainfall depth is reduced by the ratio of snow- free area to 
zonal area. At a later stage during the season, which is set to 
April 1, rain falling on snow-covered areas is allowed to 
fully contribute to runoff such that rain from entire zonal 
area is accounted. The SRM approach estimates fractional 
contributions to total discharge from (1) daily snowmelt 
from the snowpack for each elevation, (2) the melt of 
precipitation delivered as snow but to areas described as 
snow-free at the time of snowfall and (3) zonal rainfall. 
Contributions of each to total river flow are estimated prior 


to accounting for losses via runoff coefficients, introducing 
an implicit assumption that fractional contributions to runoff 
are proportional to the relative abundance of each input. 

Equation (I ) is applied to each zone of the basin, and total 
discharge is computed as sum of all zonal discharges. The 
SRM structure as applied in this study assumes a time lag 
between the daily temperature cycle and the resulting 
discharge cycle ofl 8 h. The use of this time lag is supported 
by the expected relation between lag time and basin size as 
shown in a SRM intercomparison report by the World 
Meteorological Organization (1986). Although the time lag 
parameter can be adjusted to improve synchronization of 
simulated and observed peaks of average daily flows, 
calibration of the recession coefficients provides a similar 
effect (Martinec et al., 2008). The parameters for runoff 
coefficients, lapse rate, critical temperature, recession 
coefficient and degree-day factors reflect characteristics at 
a basin scale and arc determined in this study with an MCMC 
data assimilation by maximizing the fit between the observed 
and modelled stream flows. The modelled stream flow is 
assessed for accuracy using the Nashr-Sutcliffe (NS) statistic 
and the volume difference (D v ), where volume difference 
demonstrates bias. (Supplementary information, A) 
(Martinec et al ., 2008). 

Model calibration using MCMC optimization 

A data assimilation technique known as a MCMC routine 
(Zobitz et al ., 2011) was used to explore the parameter 
spaces in the SRM. The applicability of Bayesian methods 
in hydrological studies has increased owing to their use in 
parameter estimations and uncertainly assessments (Smith 
and Marshall, 2008). The MCMC algorithm is one of 
several data assimilation techniques and uses the Metropolis 
Hastings algorithm (Metropolis et al ., 1953), which 
iteratively adjusts the model parameters to yield the best 
match between the observed and modelled strcamllow. The 
output provides a set of accepted parameter values and 
information about the posterior distributions of the 
parameters consistent with the prior information, model 
structure and observations. Additional details of the MCMC 
algorithm used for optimization are described more fully in 
Supplementary infomiation, B and in Zobitz et al. (201 1) 
and Braswell et al. (2005). 

Based on the literature and range of parameter values used 
in other snowmelt runoff models in the HKH region 
(Immerzeel et al., 2010; Butt and Bilal, 201 1; Tahir et al., 
2011), we first specified prior ranges for the parameters. 
Table II shows the parameters with the prior ranges that were 
assigned seasonally or annually across different elevation 
zones. Seasonally varying parameters were optimized for 
the snowmelt and monsoon period (June— August) and the 
snow accumulation and dry period (Sep tern ber-May ) . 
Snowmelt and rainfall runoff coefficients were varied 
seasonally and by elevations informed by Tahir et al. (201 1) 
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Table II. Initial parameter ranges used for Markov Chain Monte Carlo data assimilation 


Parameter 


Parameter ranges across each elevation zone 



Zone l 

(533- 2500 m) 

Zone 2 

(2500 4000 m) 

Zone 3 

(4000 5500 m) 

Zone 4 
(>5500 m) 

Snow runoff coefficient 

June August 

0.3 0.8 

0. 3-0.8 

0.348 

0.3-0. 8 


September May 

0.3-0. 8 

0.3 0.8 

0.3 0.8 

0,3-0. 8 

Rainfall runoff coefficient 

June August 

0.3 0.8 

0.3 0.8 

0.3 0.8 

0. 3-0.8 


Sep May 

0. 3-0.8 

0. 3-0.8 

0.3 0.8 

0,30.8 

Degree-day factor (cm/°C/day) 

Annual 

0. 3-0.6 

0. 3-0.6 

0.5 0.9 

0,5 0.9 

Lapse rale (°C/l(X)m) 

Annual 

0.5 0.7 




Critical temperature (°C) 

Annual 

0.0-= 2.0 




x coefficient 

Annual 

0.9- 1.4 




y coefficient 

Annual 

0.0 0.25 





Snowmelt and monsoonal months are June-Augus^ and Septeitlber-May is the snow accumulation and dry perk>d. Lapse rate, critical temperature and x 
and y coefficients are not varied across elevation zones. 


Lind were assigned equal prior parameter ranges. Allhough 
degree-day factors Lend to vary with different seasons, 
climate and surfaces, under similar conditions they are 
expected to increase with increasing elevation (Hock, 2003). 
Similar to previous studies (Butt and Bilal, 2011; Tahir 
et aL, 20 1 1 ), which utilize a degree-day factor that increases 
with elevation, we also increase the parameter range 
gradually across elevation zones. The elevation zone greater 
than 5500 m, which are glacier and ice covered are assumed 
to have a higher degree-day factor as the values Lend to be 
considerably higher for ice compared with snow (Hoek, 
2003). Parameters such as lapse rate, critical temperature 
Lind x and y coefficients were calibrated annually. Other 
studies using SRM have used a similar variation of 
parameters. Although most parameters can potentially vary 
across elevations and seasons, increasing parameter sets 
increase potential ranges of parameter values and thereby 
increase uncertainty (Seibert, 1997). For this study, we 
implemented the MCMC technique using three chains, each 
of which consisted of 10000 iterations. The best parameter 
values from these chains were used in the final set of 1 0 000 
iterations, which is used to determine the parameter 
distributions. The range of parameter values was deliber- 
ately set to be large at first to examine if realistic parameter 
values were selected by the MCMC technique. If the 
posterior parameter distributions were less plausible, the 
parameter ranges were logically constrained by literature 
values to be more realistic. The snowmelt and rainfall runoff 
coefficients were constrained in the ranges specified in 
Table IT based on Tahir et aL (201 1) and Butt and Bilal 
(201 1 ), degree-day factor ranges were specified on the basis 
of Hoek (2003), and prior ranges for lapse rale and critical 
temperature were specified from Tahir et aL (201 1) and Butt 
Ltnd Bilal (201 1). The SRM was run via MCMC with prior 
parameter ranges as specified in Table II for the hydrological 
years 2002 to 2006. 


Model sensitivity and uncertainty 

fn addition to model calibration, wc also carried out 
several experimental runs. The model was allowed to 
calibrate via MCMC by setting the snowmelt runoff 
coefficients to zero such that snowmelt contributions were 
neglected in the simulation of the hydrograph. Similarly, the 
rainfall runoff coefficients were also sel to zero such that 
rainfall contributions were excluded in the final simulation of 
the hydrograph. In addition to these experiments, we 
conducted additional sensitivity and uncertainty analyses 
via ensemble stream flow simulations (1000 runs) to identify 
parameters that control model performance and uncertainties 
in modelled streamllow and total input contributions. 
Parameter sensitivity and uncertainty in streamllow in the 
lime history were analyzed as 5% and 95% confidence 
intervals of the ensemble runs. First, we generated an 
ensemble of model runs ('best ensemble 7 hereafter) from the 
posterior distribution of accepted parameter sets from MCMC 
that represent the most likely set of parameters. Second, the 
sensitivity to parameters was tested by perturbing the 
distribution of a parameter by some amount and computing 
the corresponding changes in model output while preserving 
the distribution of other parameters. This is performed in this 
study by randomly drawing from the altered distribution of 
each parameter (e.g. ±30% from the mean value of a 
parameter from the posterior distribution) to generate 
ensemble urns. This allows quantifying the effects of each 
parameter or a combination of parameters on model output 
and on total input contributions to streamflow by providing 
uncertainly bounds on the output statistics. Third, an ensemble 
of parameter sets was generated by randomly sampling from 
the cumulative distribution curve constructed from the 
posterior distribution of each parameter The ensemble model 
simulation uses the full range of parameter variations and 
combinations and thereby provides overall uncertain Lies in the 
choice of model parameters generated from MCMC. 
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RESULTS 

Summary of meteorological parameters 

Long term annual temperatures (1970-2006) at the 
Taplejung station averaged 1 6.3 °C, whereas average winter 
and summer temperatures were 10 °C and 21.3 °C, 
respectively. Mann— Kendall tests indicate significant 
(p < 0.05) positive trends in both annual maximum and 
average temperatures for the Taplejung station during this 
period. Average annual precipitation at the Lungthung 
station was -2484 mm for the 1996-2006 period. On 
average, approximately 76% of the precipitation is received 
during the summer monsoon. The monthly streamflow for 
the Tamor River basin at the Majhitar station along with the 
total monthly precipitation at Lungthung station indicates a 
fingerprint of monsoonal precipitation in the Eastern 
Himalaya from June to September (Figure 3). Average 
streamflow measured at the Majhitar station averaged 
256. 2 nr/s annually during the same period, whereas the 
monsoon months had an average slreamllow of 573 m 3 /s. 
Snow cover depletion curves for the Tamor River basin 
indicate that the highest elevation zone (above 5500 m) had 
a maximum -85% snow-covered area for the 2002-2006 
period, dropping to -73% during summer (Figure 4). The 
elevation zone between 4000 and 5550 m had an annual 


maximum snow cover of -40% with summer zonal snow 
coverage falling to -12%. Elevations up to 4000 m had on 
average -9% annual maximum snow cover. 

Model parameterization 

When the prior ranges of the parameters were deliberately 
set to be large, the model calibrated well with a good fit 
between the observed and modelled streamflow; however, 
this also resulted in several models with similar performance 
and physically implausible coefficients. Logically cons- 
training the prior parameter ranges provided models with 
better fit and parameters that were physically plausible when 
varied across seasons and elevations. The posterior 
distributions of parameters after the MCMC optimization 
for all parameters for the 2002-2003 hydrological years are 
shown in Figure 5. When the model was optimized using the 
MCMC approach, some of the most sensitive parameters 
(such as x and y coefficients, critical temperature and lapse 
rate) exhibited a narrow, uni modal distribution with their 
ranges significantly reduced from the initial parameter 
ranges, The x and y coefficients required to compute the 
recession coefficient were the most sensitive parameters, as 
Lhey determined the structure of hydrograph recession and 
corresponding proportion of daily melt water that appeared 


(a) 



(c) 



(b) 



Day of year 


2002 

— 2003 

2004 

2005 

* 2006 


Figure 4. Snow cover depletion curves from the Moderate Resolution Imaging SpecLmradiometer 8-day maximum snow cover 500-m resolution product in the 
Tamor River basin from 2002 to 2006 in (a) /one 2, elevation range Trom 250CMO00 m, (b) zone 3, elevation range from 4000-5500 m and (c) zone 4, elevation 

range greater Lhan 5500 m, 
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Figure 5, Posterior parameter distributions tor the 2002-2003 hydrological years as histograms 
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immediately in the runoff. This shows that the MCMC 
approach performs well in constraining the most sensitive 
parameters of the model. Once a narrow posterior 
distribution is established for the sensitive parameters 
during the MCMC approach, changing other parameters 
does not result in a significant increase in model 
performance. There tore, some of the remaining parameters 
such as snowmelt and rainfall runoff coefficients exhibited 
broad posterior distributions, and large standard deviations 
where good simulations of the observed streamflow were 
obtained over broad ranges (Figure 5 and Table 111). The 
posterior distributions of runoff coefficients for both rainfall 
and snowmelt are comparable across types of precipitation 
and all simulated years. These coefficients represent the 
fraction of liquid water, from either snowmelt or rainfall, 


that contributes to runoff and are expected not to vary as 
much across types of precipitation. Some of the observed 
differences in these coefficients within a single zone could 
be a result of the trade-off among parameters after the most 
sensitive parameters have been constrained well. Inadequate 
precipitation data from higher elevations zones in Lhe 
mountainous regions further hindered our ability to 
constrain the runoff coefficients. Parameter estimates tend 
to be stable across years, though Tcrit showed large 
variation, but also had a large spread. Degree-day factor 
estimates tended to be the highest for elevations greater than 
4000 m (zones 3 and 4), which had the highest fraction of 
snow-covered area. The peaks in posterior distributions for 
the degree-day factor correspondingly rise from zone 1 to 
zone 3 elevation ranges. 
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Table HI. Mean and standard deviation of parameters from the posterior distribution 


Parameter 



Mean and standard deviation of parameter across zones 


Year 

Zone 1 

(533- 2500 m) 

Zone 2 

(25(X) 4fXX) m) 

Zone 3 

(4000 5500 m) 

Zone 4 
(>5500 m) 

Cs 

June August 

2002 

0.55 ±0.15 

0.59 ±0.1 3 

0.72 ± 0.07 

0.56 ±04 3 



2<M)3 

0.46 ±0.1 3 

0.44 ±0.11 

0.47 ±0.1 2 

0.52 ±045 



2004 

0.54 ±0.14 

0.54 ±0.14 

0.68 ±0.11 

0.55 ±044 



2005 

0.55 ±0.14 

0.55 ±0.1 5 

0.69 ±0.11 

0.55 ±045 


September May 

2002 

0.58 ±0.14 

0.61 ±0.11 

0.59 ±0.1 3 

0.54 ±045 



2003 

0.60 ±0.14 

0.56 ±0.1 3 

0.71 ±0.07 

0.56 ±04 5 



2004 

0.70 ±0.09 

0.71 ±0.07 

0.67 ±0.11 

0.54 ±04 4 



2005 

0.59 ±0.14 

0.43 ±0.11 

0.40 ±0.0/8 

0.54 ±04 4 

Cr 

June August 

2002 

0.73 ±0.05 

0.74 ±0.05 

0.69 ±0.09 

0.55 ±045 



2003 

0.67 ±0.11 

0.56 ±0.13 

0.54 ±0.1 5 

0.54 ±04 4 



2004 

0.71 ±0.05 

0.68 ±0.09 

0.59 ±0.14 

0.55 ±04 4 



2005 

0.76 ±0.04 

0.74 ±0.08 

0.76 ±0.05 

0.56 ±044 


September May 

2002 

0.69 ±0.07 

0.68 ±0.09 

0.65 ±0.10 

0.55 ±044 



2003 

0.37 ± 0.07 

0.43 ±0.09 

0.40 ±0.08 

0.54 ±04 4 



2004 

0.71 ±0.07 

0.70 ±0.08 

0.54 ±0.1 4 

0.56 ±04 4 



2005 

0.69 ±0.07 

0.70 ±0.08 

0.72 ±0.07 

0.51 ±043 

a 

Annual 

2002 

0.48 ±0.08 

0.50 ±0.0 

0.61 ±0.08 

0.68 ±0402 



2003 

0.42 ±0.08 

041 ±0.08 

0.79 ±0.0 

0.69 ±041 



2004 

0.53 ±0.05 

0.53 ±0.06 

0.81 ±0.05 

0.70 ±022 



2005 

0.45 ±0.09 

0.44 ±0.09 

0.80 ±0.08 

0.71 ±042 

LR 

Annual 

2002 

0.64 ±0.01 






2003 

0.51 ±0.01 






2004 

0.68 ±0.01 






2005 

0.54 ±0.01 




Tent 

Annual 

2002 

1.26 ±0.28 






2003 

1.1 6± 0.49 






2004 

1.83 ±0.13 






2005 

0.94 ±0.62 




x coefficient 

Annual 

2002 

1.34 ±0.03 






2003 

1.26 ±0.01 






2004 

1.29 ±0.02 






2005 

1.11 ±0.01 




y coefficient 

Annual 

2002 

0.08 ±0.01 






2003 

0.06 ±0.01 






2004 

0.07 ±0.01 






2005 

0.03 ±0.01 





Cs, snowmelt runoff coefficient; Cr, min fall rnnoff C Deficient; a, degree day factor; LR, lapse rale. 


The SRM is usually calibrated manually where known 
ranges of parameter values are maintained while eh tinging 
other parameters by trial and error. This type of calibration 
can be inefficient in exploring the parameter space and 
resulted in an average NS statistic of 0.64 and volume 
difference of 10% across all simulated years for this study. 
Model calibration using MCMC optimization had a greater 
overall improvement with NS statistic of -0.84 averaged 
across all simulated years. Figure 6 (a, e, c and g) shows the 
modelled stream flow against observed slreamflow using the 
MCMC approach for four hydrological years (utilizing 
observed precipitation at the Lungthung station). The SRM 
is able to capture the daily slreamflow in the annual 
hydrology of the basin; however, the model undcres Li males 
high discharge events, particularly during the monsoonal 


months. Across all simulated years, the observed and 
modelled slreamflow show a linear association with some 
bias during high discharge periods [Figure 6 (b, d, f and h)]. 
The model simulated streamflows have a higher NS statistic 
in the years where the summer peak discharge is lower 
relative to other years (Table TV and Figure 6). 

Model sensitivity and uncertainty 

The results of the experiment runs indicate that with the 
snowmelt runoiT coefficients (Cs) set to zero, the model 
optimized as well as when Cs values were also included 
through the MCMC technique. The rainfall runoff 
coefficients (Cr) are parameterized to higher values, and 
Lhe model is able to simulate the hydrograph by optimizing x 
and y coefficients with model performance equivalent to that 
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Figure 6. Time scries and scatter plots of observed and modelled streamflow using observed precipitation at Lungthung station for 2002 2003 (a, b), 

2003 2004 (c, d), 2004 2005 (c, f) and 2005 2006 (g, h) 
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with the best estimates of parameters as in Table IV. The first 
ensemble parameter sets (best ensemble) randomly selected 
from the posterior distribution output from MCMC did not 
exhibit any significant variations in model fit and perfor- 
mance. This is to be expected as the parameter sets are 


accepted by the MCMC approach. Because the most 
sensitive parameters have narrow, unimodal distributions, 
similar model fit and performance is achieved within a broad 
distribution of other parameters. However, when the 
distribution of the most sensitive parameters was altered 
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Table IV. Evaluation of snowmelt runoff model performance wilh Markov Chain Monle Carlo technique and using observed 
precipitation at Lungthung station and gridded Asian Precipitation Highly Resolved Observational Data Integration towards 

Evaluation precipitation data 

Hydrological Model efficiency (wilh observed precipitation) Model efficiency (with APHRODITE precipitation) 

Year 



Nash -Sutcliffe 
metric 

Coefficient of 
determi nation 

Difference of 
volume (D v %) 

Nash Sutcliffe 
metric 

Coefficient of 
detenu inalion 

Difference 
of volume (D v %) 

2002 2003 

0.84 

0.84 

1.9 

0.73 

0.79 

13.1 

2003 2004 

0.85 

0.85 

5.2 

0.82 

0.87 

11.9 

2004 2005 

0.80 

0.80 

2.7 

0.78 

0.80 

5.10 

2005 2006 

0.87 

0.87 

0.7 

0.85 

0.85 

3.12 


APHRODITE* Asian Precipitation — Highly Resolved Observational Data Integration towards Evaluation. 


slightly from their posterior distributions in the second 
ensemble parameter sets, there were notable changes in the 
model fit of the ensemble tuns. Changing the x parameter 
randomly within ±5% of its mean value of its posterior 
distribution resulted in sizeable changes in model fit and 
performance (NS ranges from 0.3-0. 8 and D v ranges from 
□6% to 12%). Similarly, changing the y parameter also 
resulted in notable changes in model fit and performance 
(NS ranges from 0.7-0. 8 and D v ranges from 8% to 7%). 

When the parameter distributions ate altered randomly 
within ±10% of mean value of the posterior distribution, 
degree-day factor and critical temperature did not affect the 
model fit and total input contributions relative to the best 
ensemble. When lapse rate w r as varied within ± 1 0% of mean 
value of its distribution, the average total input contributions 
of snowmelt were 25.8 ±4.4%, 4.2 ± 1 . 9 % for precipitation 
as snow on previously snow-free areas, and 70.0 ±5.8%: for 
contributing precipitation across the 2002—2006 period. 
Figure 7a shows a larger 5th and 95th percentile envelope 
for simulated streamflows using lapse rate, indicative of the 
model's sensitivity to this parameter relative to the narrower 
envelope for critical temperature in Figure 7b. 


The final ensemble runs where the parameters are 
randomly sampled from the cumulative distribution curve 
capture the uncertainly over the full range of parameter 
variations from the posterior distribution. Figure 8 shows the 
uncertainty bounds from the ensemble streamilow, which 
are large enough during the spring snowmelt period to 
include most of the observations but not large enough 
particularly during the summer monsoon months. The final 
ensemble streamflow simulations were also used to generate 
average total input contributions of rainfall and snowmelt 
for all simulated years and standard errors of these 
contributions (Figure 9). The average proportion of 
snowmelt in the runoff in the Tamor River basin for the 
2002-2006 period is 29.7 ±2.9%, which includes 
4.2 ± 0.9% from new snow onto previously snow-free areas, 
whereas 703 ± 2.6% is attributed to rainfall contributions. 
The cumulative runoff indicates that most of the snowmelt 
contribution occurs between April and August, Although 
there is interannual variability in the contribution of 
snowmelt runoff, elevation zone 3 (4000^5500 m) contrib- 
utes on average 56.9 + 3.6% of all snowmelt input and 
28.9+ I J% of all rainfall contribution to runoff. Elevation 


(a) (b) 



Month Month 


Figure 7, Slreamtluw ensemble simulations (2002-2003) by randomly varying the distribution oT the parameters by ±10% of its mean posterior 
distribution. Black circles are observations and the dark cyan lines depict the 5lh and 95 th percentiles of the ensemble. Where lines for the 5lh and 95th 

percentiles overlap only a single line is visible 
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Figure 8. Slroamflow ensemble simulations using observed precipitation at Lunglhung station for Tamor River basin Tor hydrologic years (a) 2002-2003, (b) 
2003-2004, (c) 2004—2005 and (d) 2005-2006. Black circles are observations and the cyan lines depict the 5lh and 95lh percentiles from the ensemble 
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Figure 9. Cumulative curves of computed daily snowmelt depth, melted precipitation in the form of snow and rainfall depths for the (a) 2002-2003, 
(h) 2003—2004, (c) 2004-2005 and (d) 200S-2006 hydrological years for the Tamor basin 
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zones 1 (5 33-2500 m) and 2 (2500-4000 m) contribute calibrated using MCMC captured the observed streamflow 

on average 67.9 ±2.7% of Lfie all annual rainfall input to fairly well (Figure 10). Modelled stream flow using observed 

the basin. station precipitation is also plotted alongside, which 

performs better on the basis of evaluation metrics and 
Model results using APHRODITE precipitation simulation of peak discharges (Table IV). Although the NS 

The modelled stream flow obtained using APHRODITE coefficient value was -0.80 on average across all 

precipitation data over the hydrological periods and simulations, the model failed to capture peak discharge 





Month 



Month 





Figure 10. Time series and scalterploLs of observed and modelled stream!] ow using Asian Precipitation - Highly Resolved Observational Data Integration 
towards Evaluation precipitation (black) and observed station precipitation (blue) for 2002-2003 (a, b), 2003-2004 (c, d), 2004—2005 (e, f) and 2005-2006 (g, h) 
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periods contributing to difference in volume (D v ) of on 
average 8.3%. Generally, the modelled hydrology using 
basin-averaged APHRODITE precipitation failed to capture 
high discharge periods across all years when compared to 
the simulations using gauged precipitation as input data. 
Overall, use of the APHRODITE precipitation data yielded 
similar model skill in simulating the hydrograph but 
estimated a smaller total contribution from precipitation 
input to the basin (more from snowmelt) when compared 
with model simulations using the observed station precip- 
itation data. 


DISCUSSION 

This study provides an application and evaluation of the 
SRM in a mountainous basin in the eastern Nepalese 
Himalaya. In this study, the SRM was coupled to a MCMC 
data assimilation approach to examine associated param- 
eter sensitivities, uncertainties and model performance in 
the Tamor River basin for the 2002-2006 period. A 
majority of the studies involving snowmelt runoff models 
use parameter values from existing studies, and therefore, 
parameter calibration, types of input data and equifinality in 
parameter sets have not been fully explored (particularly for 
the Himalayan region). Additionally, lack of observational 
data is problematic for transfer of parameters across space 
and time and may introduce added uncertainly. Uncertainty 
in a model is usually a joint outcome of uncertainties in 
model structure, input data errors and parameter errors, 
which leads to difficulty in approximating the behaviour of a 
model consistently (Vrugt et al ., 2005). Although any 
parameter optimization technique can only reduce uncer- 
tainty in the parameters, optimization and data assimilation 
methods such as the MCMC may be better for accounting 
important sources of uncertainty. 

Overall, parameter optimization using MCMC increased 
the fit between the observed and modelled stream flows and 
further reduced the volume difference error. This study also 
showed that prior ranges of parameter values are important 
in constraining model performance and decreasing param- 
eter uncertainty. When the parameter ranges were logically 
constrained or infonned from estimates from previous 
studies and then optimized using MCMC, the output 
parameters value ranges were more relevant and physically 
plausible. Model evaluation through additional experiment 
runs indicate that despite overall good performance of SRM, 
there might be lack of confidence in model structure in 
simulating physical realism of some of the processes. For 
instance, the model was able to reproduce the hydrograph 
well even when snowmelt runoff contributions were 
excluded from simulations. As the contributing rainfall to 
the Tamor basin is high, relative to the snow component, 
model parameters were able to compensate for the absence 


of snowmelt contributions and simulate the annual 
hydrograph with relatively high fit and performance. This 
lest identifies an important model structural uncertainly that 
is not currently reflected in the reported uncertainty on 
forecasts and the associated snowmelt contribution. It also 
identifies that the hydrograph does not constrain estimates of 
the fractional contributions of total outflow coming from 
snowmelt versus rainfall. Instead, the model 1 s estimates of 
the snowmelt contribution derive mainly from volume 
considerations based on the snow-covered area, the degree- 
day model of melt and rainfall depths. This lends important 
new insight into the practical outcomes of SRM application 
that may not have been given full attention in past 
applications for similar purposes. This partly reflects a 
limitation of a conceptual lumped model such as the SRM. 
However, there are some clear advantages to using only a 
few parameters for modelling, particularly when there is 
such a lack of observations and strong seasonal variation in 
the nature of parameters across a topographically challeng- 
ing environment. When more complex physically based 
hydrological models are applied to such environments, 
calibration can be further challenged owing to the large 
number of parameters involved. Furthermore, it is not clear 
that such an improvement in mechanistic realism would 
translate into improved model skill in matching the 
hydrograph or in representing volume contributions to total 
basin surface water out flow. 

The posterior distribution output from MCMC indicates 
that the coupled approach was successful in constraining the 
ranges of sensitive parameters in the model. Model 
sensitivity using ensemble runs indicated that the x 
coefficient, y coefficient and lapse rate are the most sensitive 
parameters in the model. The recession coefficient (&), which 
is a function of x and y coefficients, is important as it 
represents an inertia or memory in stream flow and also 
controls the amount of the daily melt water that appears in the 
mnoff, whereas lapse rate is important in lapsing the zonal 
temperature to compute daily zonal snowmelt depth. The 
narrow, uni modal distributions of these parameters provide 
an indication of the degree to which these parameters were 
constrained by the MCMC approach. Parameters Lhal are 
well constrained tend to have small standard deviations with 
well-defined unimodal distributions, whereas poorly 
constrained parameters exhibit flat distributions with large 
standard deviations (Braswell et al. f 2005). Once the 
sensitive parameters are constrained to a narrow distribution, 
other parameters that are not as sensitive provide a similar 
model fit and performance resulting in broad posterior 
ranges. A study by Seibert (1997) also found that only few of 
the parameters were well defined when the Hydro logi ska 
By ran s Vattenbalansavdelning model, a conceptual degree- 
day model, was parameterized over central Sweden. This is 
also suggestive of trade-offs among parameters leading to 
broader posterior distributions where different parameter 
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sets may be acceptable in reproducing the observed 
hydrology within a basin. 

The overall model parameter uncertainty was character- 
ized as standard errors in the total input contributions and 
uncertainty in time history of streamflow as 5% and 95% 
confidence intervals (Figure 8). Because the uncertainty 
bounds are not large enough to include all of the 
observations, parameter uncertainly alone cannot explain 
the total error in the model. If the bounds were large enough, 
parameter variability alone could account for overall 
uncertainty by compensating other sources of error intro- 
duced by model structure and input data (Bias one et al . , 
2008). As most of the observations outside of the 
uncertainty bounds are during the summer months, the 
unaccounted uncertainty could arise from uncertainty in 
input precipitation data. 

On the basis of the evaluation metrics, use of APHRO- 
DITE gridded precipitation product provided comparable 
skill relative to the discharge modelled using observed 
station observations. APHRODITE data have been shown 
to be effective in simulating and forecasting discharges 
using a snowmelt model in other high-altitude ranges such 
as the Karakoram region (Tahir et al., 2011). The results 
indicate that the total rainfall contributions to runoff is lower 
when APHRODITE is used for calibrating the model. This 
is consistent with past studies that indicated underestimated 
total annual precipitation and high discharges when 
APHRODTTE was used in snowmelt runoff modelling in 
the Hunza River basin (Tahir et al. y 2011). Nonetheless, the 
findings from our study using APHRODITE indicate that 
these data can be important inputs for precipitation in basins 
lacking ground observations. We also find, however, that 
uncertainties in the available input data such as precipitation 
are equally as important as uncertainty introduced by model 
parameters. Model parameters calibrated using an input 
precipitation that either underestimates or overestimates 
actual precipitation will compensate for such error when 
matching the hydrograph (Pcllieciotli et al . , 2012) but may 
still alter conclusions about fractional contributions from 
snowmelt and rain. 

This study shows a total snowmelt contribution to be 
29.7 ±2.9% of annual discharge averaged across the 
2002—2006 hydrological years for the Tamor basin. 
Bookhagen and Burbank (2010) showed that Eastern 
Himalayan catchments draining into the Bay of Bengal 
received -80% of their annual rainfall during monsoon and 
less than 20% from snowmelt. For comparison, snowmelt 
contributions to annual runoff in the Satluj basin 
(22275 km 2 ) in the western Himalayan region has been 
estimated to be up to -59%, with the summer contributions 
averaging -75% for the late 1980s period (Singh and Iain, 
2003). Another study modelling snow and glacial melt in the 
Dudh Kosi River basin (3712 km 2 ) estimated total melt 
contributions to be 34% (Nepal etal . , 2013). This study also 


found the 4000-5500 m elevation range to be the most 
important with regard to its snowmelt contributions to 
annual runoff. Biggs and Whitaker (2012) also introduced 
the concept of critical zone using the example of Merced 
River basin, California, which can be important in 
anticipating the impacts of regional climate change on the 
Liming and volume of runoff events. 

A logical incremental improvement on this work would 
be to incorporate the present and future role of glaciers in 
the runoff. Previous studies in the Himalayan region (e.g. 
Tahir et al . , 2011) have used the SRM to include glacier 
melt by calibrating the model using a higher degree-day 
factor for the higher elevations. Although this SRM study 
does not explicitly incorporate glaciated regions, they are 
partly taken into account through the use of snow depletion 
curves. The decline of snow cover depletion curves levels 
out during the later portion of the snowmelt season, thereby 
indicating the presence of permanent snow or glaciers (e.g. 
zone 4). Model modifications would involve incorporation 
of glacier-covered areas at the highest elevations and 
calibration of model for parameters such as degree-day 
factor that are different for snow and ice. Additional input 
data such as estimates of snow cover data, calibration 
against multiple variables such as inclusion of glacier mass 
balances or fixing known parameters from field observations 
are usually the options to further narrow the posterior 
distributions and better constrain these parameters that vary 
seasonally across elevations (Boudhar et al ., 2009; 
Pclliceiotti et al., 2012). Given the limited availability of 
data and the parsimoniousness of this snowmelt runoff 
model, obtaining values through held observations or fixing 
parameters based on values reported from previous short 
field campaigns are potential options to improve the model. 


SUMMARY 

This study has shown that a snowmelt runoff model based 
on degree-day factors has skill in simulating daily 
streamflow in a mountainous environment with limited 
coverage of hy drome Leorological measurements. Coupling 
the SRM with an MCMC method provided an improved fit 
to observed streamflow s but still failed to capture peak 
discharge during the summer monsoon months. This is 
likely due in part to the lack of precipitation data coverage 
throughout the basin but could also result from model 
structural errors. Model performance in simulating the 
hydrograph is strongly sensitive to recession coefficients 
and the lapse rate, but exhibited little sensitivity to runoff 
coefficients, the critical temperature that determines the 
phase of precipitation (snow or rain) and the degree-day 
factor that estimates melt depth from degree days. 
Correspondingly, the MCMC approach to parameter 
determination yielded narrow distributions only for those 
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parameters to which the model -estimated hydro graph is 
sensitive. The study provides a useful guide for how to 
constrain parameters, provide uncertainty bounds in snow- 
melt contributions to runoff, analyze effects of input 
precipitation and examine model uncertainty in Himalayan 
basins. The Himalayan terrain, meteorological conditions 
and lack of ground-based observations continue to challenge 
snowmelt runoff and hydrologic modelling in the region. By 
testing results with two alternative precipitation datasets, we 
noted how both the available input data and the parameter 
sets introduce uncertainties, and methods to thoroughly 
characterize and minimize these uncertainties is important. 
The need for further hydrological research on snow and 
glacial melt and, moreover, strengthening of hydrometeo- 
rological observation network for the HKH region should 
not be overlooked. Overall, this study shows that a 
combination of known parameter ranges along with an 
iterative model-data fusion technique such as the MCMC 
can provide a calibration and optimization method to 
explore the parameter space efficiently, obtain reliable 
parameter estimates and quantify uncertainty in model 
results and associated inferences. 
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